1  Aires protégées

1.1 Données de l’association Vahatra

Les études sur les aires protégées s’appuient fréquemment sur la base WDPA (World Database on Protected Area), consultable en ligne sur https://protectedplanet.net. On s’aperçoit dans le cas de Madagascar que cette base de données comporte de nombreuses erreurs (qu’on étudiera plus bas). La base rassemblée par l’association Vahatra dans le cadre de la monographie qu’elle a coordonnée sur l’ensemble des aires protégées terrestres malgaches semble beaucoup plus fiable (Goodman et al. 2018). Les données en question sont disponibles sur le portail https://protectedareas.mg avec une licence creative commons (CC-BY).

Le bloc de code ci-dessous (cliquer sur “code” pour visualiser), présente la séquence d’opérations réalisées pour préparer les données.Pour comprendre certaines opérations contenues dans le bloc de code, il est utile d’être familier de la syntaxe de R et des packages du tidyverse. Voir le chapitre Chapter 10.

Code

library(tidyverse)
library(lubridate)
library(sf)
library(tmap)
library(geodata)
library(cowplot)
library(wdpar)
library(gt) # Pour faciliter le rendu des tableaux (et ils sont jolis)

# Le shapefile est composé d'une série de fichiers, (.shp, .dbf, .prj, .shx)
# qui doivent avoir le même nom et être au même endroits pour être ouverts en
# même temps. Comme souvent, ils sont compressés ensemble dans un fichier zip.
# On commence par dézipper (décompresser) ce fichier.
unzip("data/Vahatra98AP.zip", exdir = "data/Vahatra")
# On importe dans R en pointant vers le fichier .shp, mais c'est bien toute la
# collection de fichiers homonymes .shp, .dbf, .shx qui est chargée.
AP_Vahatra <- st_read("data/Vahatra/Vahatra98AP.shp", quiet = TRUE) %>%
  # Il manque la projection (pas de fichier .prj), on la spécifie à la main
  st_set_crs("EPSG:4326") # EPSG 4325 = WSG 84 = le standard pour le web

# L'option ci-dessous est un peu cryptique : des caractéristiques topologiques
# de la carte source sont incompatibles avec la possibilité d'avoir des objets
# sphériques dans sf. Cela disparait si on désactive cette possibilité
sf_use_s2(FALSE) 

# Identification des dates ----------------------------------------------------
# Cette section est un brin complexe, à base de manipulation de chaînes de 
# caractères et de dates

# Détecte les dates écrites 2 avril 2020 ou 02 avril 2020, etc.
date_ecrite <- "[:digit:]{1,2} [:alpha:]* [:digit:]{4}"
# Détecte les dates écrites 02/04/20 ou 02.04.20 ou 02.04.2020, etc.
date_abrev <- "[:digit:]{2}[:punct:][:digit:]{2}[:punct:][:digit:]{2,4}"
# Des années écrites à 2 chiffres
date_ecrite_an_abrev <- "[:digit:]{1,2} [:alpha:]* [:digit:]{4}"
# Détecte l'une ou l'autre des formes précédentes
toute_date <- paste(date_ecrite, date_abrev, date_ecrite_an_abrev, sep = "|")
# Détecte une mention d'année seule : 1984, 2015, etc.
annee_seule <- "[:digit:]{4}"
# Détecte les formes indicatrices d'un changement
mention_changement <- "Changement|changement|anciennement|actuel|auparavant"
# Une fonction qui traduit les dates écrites en toutes lettre du français à 
# l'anglais (pour les parser ensuite car ça ne fonctionne qu'en anglais)
trad_dates <- function(date_fr) {
  str_replace_all(date_fr,
                  c("janvier" = "January",
                    "fevrier" = "February",
                    "mars" = "March",
                    "avril" = "April",
                    "mai" = "May",
                    "juin" = "June",
                    "juillet" = "July",
                    "aout" = "August",
                    "septembre|setembre" = "September",
                    "octobre" = "October",
                    "novembre" = "November",
                    "decembre|decmbre" = "December"))
}
# Cette fonction remplace 01.04.58 par 01.04.1958 et marche avec . ou /
# On indique avec limite le nombre d'année où on considère que c'est 1900 vs 2000
complete_annee <- function(date_abrev, limite = 20) {
  if (str_detect(date_abrev, "([:punct:])([:digit:]{2})[:punct:]?$")) {
    date_abrev <- str_remove(date_abrev, ":punct:]?$")
    if (as.numeric(str_extract(date_abrev, "[:digit:]{2}$")) > limite) {
      date_abrev <- str_replace(date_abrev, 
                                "([:punct:])([:digit:]{2})[:punct:]?$", "\\119\\2")
    } else {
      date_abrev <- str_replace(date_abrev, 
                                "([:punct:])([:digit:]{2})[:punct:]?$", "\\120\\2")
    }
  }
  return(date_abrev)
}
# La fonction précédente est unitaire, on la transforme pour qu'elle s'applique à une liste.
complete_liste_dates <- function(liste_dates) {
  map(liste_dates, complete_annee)
}

AP_Vahatra <- AP_Vahatra %>%
  # On extrait les dates des champs de texte
  mutate(date_creation = str_extract_all(creation, toute_date), 
         # Une date a un format incohérent, on la recode à la main
         date_creation = ifelse(creation == "Créée le 07 aout 04",
                                "07 aout 2004", date_creation),
         date_creationA = map(date_creation, 1), # La 1ère date
         date_creationB = map(date_creation, 2)) %>% # Si 2 dates, la seconde
  # On traduit les mois en anglais pour une conversion au format date
  mutate(across(c("date_creationA", "date_creationB"), trad_dates)) %>%
  mutate(across(c("date_creationA", "date_creationB"), complete_liste_dates)) %>%
  mutate(across(c("date_creationA", "date_creationB"), dmy)) %>%
  mutate(date_creation = case_when(is.na(date_creationB) ~ date_creationA,
                                    date_creationA > date_creationB ~ date_creationB,
                                    date_creationA <= date_creationB ~ date_creationA),
         date_modification = case_when(is.na(date_creationB) ~ date_creationB,
                                       date_creationA < date_creationB ~ date_creationB,
                                       date_creationA >= date_creationB ~ date_creationA),
         # On repère si il y a eu un changement de statut ou de frontières
         mention_changement = str_detect(creation, mention_changement)) %>%
    # On enlève les colonnes inutiles
  select(-date_creationA, -date_creationB) %>%
  # On place les colonnes créées à gauche pour les inspecter facilement
  relocate(date_creation:mention_changement, .after = creation) 

# Après une vérification manuelle, on remarque les données de l'association Vahatra comportent des mentions incomplètes pour certaines aires, qui n'ont pas été extraites:


# Lokobe : 31 décembre 1927
AP_Vahatra <- AP_Vahatra %>%
  mutate(date_creation = case_when(nom == "Lokobe" ~ ymd("1927-12-31"),
                                   nom == "Mantadia" ~ ymd("1989-01-11"),
                                   TRUE ~ date_creation),
         date_modification = case_when(nom == "Bemaraha" ~ ymd("2011-07-06"),
                                       nom == "Lokobe" ~ ymd("2011-07-06"),
                                       nom == "Tsaratanana" ~ ymd("2011-07-06"),
                                       nom == "Pic d'Ivohibe" ~ ymd("2015-04-28"),
                                       nom == "Mantadia" ~ ymd("2002-08-07"),
                                       TRUE ~ date_modification)) %>%
  st_make_valid() # fiabilise qu'il n'y a pas d'erreurs topologiques
# dir.create("AP_Vahatra")
# st_write(AP_Vahatra, "out/AP_Vahatra.shp")
# writexl::write_xlsx(st_drop_geometry(AP_Vahatra), "AP_Vahatra.xlsx")

Le bloc de code suivant génère une carte interactive. On a également inclus des lignes de code qui permettent de formater la carte joliment pour un rendu figé (pdf/LaTeX, html statique, word), mais ce code est “commenté”, c’est-à-dire qu’on a placé des dièses au début de chaque ligne, de sorte qu’il ne s’exécute pas (R n’exécute jamais ce qui se trouve à droite d’un # sur une ligne). Pour plus de détails sur la manière dont on produit des cartes, voire l’annexe : Cartes simples en R

Code

if (file.exists("data/contour_mada.rds")) {
  load("data/contour_mada.rds")
} else {
  contour_mada <- gadm(country = "Madagascar", resolution = 1, level = 0,
                     path = "data/GADM") %>%
  st_as_sf()
# On enregistre contour_mada pour s'en servir par la suite
save(contour_mada, file = "data/contour_mada.rds")
}

# On génère un rendu cartographique
tmap_mode("view") # En mode interactif
tm_shape(contour_mada) +
  tm_borders() +
  tm_shape(AP_Vahatra) + 
  tm_polygons(col = "cat__iucn", alpha = 0.6, title = "Catégorie IUCN") +
  tmap_options(check.and.fix = TRUE) # +

Code

# Les dièses en début de ligne font que ce qui suit ne s'exécute pas.
# La suite est uniquement pour les rendus fixes (tmap_mode = "plot"), p. ex. pour les pdf
  # # NB : on note les positions en majuscules quand on veut coller aux marges
  # tm_credits("Sources: WDPA et GADM", position = c("RIGHT", "BOTTOM"),
  #            size = 0.6) +
  # tm_layout(main.title = "Aires protégées de Madagascar",
  #           # NB : position en minuscules pour laisser un espace avec la marge
  #           main.title.position = c("center", "top"),
  #           main.title.size = taille_titres_cartes,
  #           legend.position = c("left", "top"),
  #           legend.outside = TRUE)

On peut également réaliser un graphique qui présente l’historique de création des aires protégées. Pour plus de précisions sur la manière de produire des graphiques en R, voir l’annexe correspondante.

Code

# On ordonne les nom d'aires protégées dans l'ordre de leur séquence de création
ordre_chrono_AP <- AP_Vahatra %>%
  arrange(desc(date_creation), desc(nom)) %>%
  pull(nom)
# On transforme le champ "nom" de caractère, à une catégorisation ordonnée où
# l'ordre correspond 
AP_Vahatra_carte <- AP_Vahatra %>%
  mutate(nom = factor(nom, levels = ordre_chrono_AP),
         cat_taille = case_when(hectares > 300000 ~ 2,
                                hectares > 150000 ~ 1.5,
                                hectares >  50000 ~ 1,
                                             TRUE ~ 0.5)) %>%
  rename(`Catégorie IUCN` = cat__iucn)

# On crée un graph pour les anciennetés
graph_gauche <- AP_Vahatra_carte %>%
  ggplot(aes(x = date_creation, xend = ymd("2022-10-01"), y = nom, yend = nom, 
                   color = `Catégorie IUCN`)) +
  geom_segment(size = 2) +
  ggtitle("Ancienneté") +
  theme(axis.ticks.y = element_blank(),
        legend.position = "none",
        axis.title = element_blank(),
        axis.text.x.bottom =  element_text(angle = 45, hjust = 1),
        axis.text.x.top = element_text(angle = -45, hjust = 1)) + 
  scale_x_date(sec.axis = dup_axis())

graph_droite <- AP_Vahatra_carte %>%
  ggplot(aes(x = 0, xend = hectares/100, y = nom, yend = nom, 
                   color = `Catégorie IUCN`)) +
  geom_segment(size = 2) + 
  ggtitle("Surface (km2)") +
  theme(axis.text.y = element_blank(),
        axis.title = element_blank(),
        axis.text.x.bottom =  element_text(angle = 45, hjust = 1),
        axis.text.x.top = element_text(angle = -45, hjust = 0),
        legend.position = "none") + 
  scale_x_continuous(sec.axis = dup_axis())

legende <- get_legend(graph_gauche  +
                        guides(color = guide_legend(nrow = 1)) +
                        theme(legend.position = "bottom"))

# On colle les deux
graphs <- plot_grid(graph_gauche, graph_droite, rel_widths = c(2.2, 1),
          nrow = 1)
plot_grid(graphs, legende, ncol = 1,
          rel_heights = c(1,.1))

Il faut aussi s’assurer qu’on filtre bien les entités analysées selon un critère pertinent. Actuellement, on exclut les aires marines. Il pourrait toutefois sembler utile d’écarter les aires dont le statut de protection est considéré comme trop faible. Il pourrait aussi être pertinent de ne garder que les aires protégées comportant un niveau minimum de couvert forestier : autrement, cela signifie que la forêt n’est pas un habitat pertinent pour les écosystèmes que la démarche de conservation cherche à protéger dans cette aire.

1.2 World Database on Protected Areas

On commence par télécharger et présenter ces données.

Code

# Il y a un bug dans l'environnement Quarto, on installe la version de dévelop-
# pemment car ça résoud le problème.
# Downloading Protected areas from Madagascar
WDPA_Mada <- wdpa_fetch("Madagascar", wait = TRUE,
                      download_dir = "data/WDPA") 

# TODO: inclure graph et description des données.

1.3 Comparaison des données Vahatra et WDPA

On commence par visualiser les différences spatiales entre les polygones, en affichant les 10 qui sont les plus différents entre les WDPA et Vahatra.

Code

# On harmonise les noms qui sont parfois notés différemment entre les sources
AP_Vahatra <- AP_Vahatra %>%
  mutate(nom_wdpa = case_when(
    nom == "Corridor Forestier Bongolava" ~ "Corridor forestier Bongolava",
    nom == "Ranobe PK32" ~ "Ranobe PK 32",
    str_detect(nom, "Ambositra-Vondrozo") ~ "Corridor Forestier Ambositra Vondrozo",
    nom == "Réserve deTampolo" ~ "Réserve de Tampolo",
    nom == "Bombetoka Beloboka" ~ "Bombetoka Belemboka",
    nom == "Ampananganandehibe-Behasina" ~ "Ampanganandehibe-Behasina",
    nom == "Forêt Sacrée Alandraza Analavelo" ~ "Analavelona", # vérfié sur carte : les mêmes
    nom == "Réserve speciale Pointe à Larrée" ~ "Réserve spéciale Pointe à Larrée", 
    nom == "Vohidava-Betsimalaho" ~ "Vohidava Betsimalao", 
    nom == "Anjanaharibe Sud" ~ "Anjanaharibe_sud",
    nom == "Iles Radama/Sahamalaza" ~ "Sahamalaza Iles Radama",
    nom == "Kalambatritra" ~ "Kalambatrika",
    nom == "Mananara-Nord" ~ "Mananara Nord",
    nom == "Kirindy - Mitea" ~ "Kirindy Mite",
    nom == "Midongy du Sud" ~ "Befotaka Midongy", # Vérifié sur la carte
    nom == "Montagne d'Ambre/Forêt d'Ambre" ~ "Montagne d'Ambre",
    nom == "Tsimanampesotsa" ~ "Tsimanampesotse",
    nom == "Pic d'Ivohibe" ~ "Ivohibe",
    nom == "Forêt Naturelle de Petriky" ~ "Forêt Naturel de Petriky",
    nom == "Tsingy de Namoroka" ~ "Namoroka",
    nom == "Réserve de Ressources Naturelle Mahimborondro" ~ "Mahimborondro",
    str_detect(nom, "Complexe Tsimembo Manambolomaty") ~ "Complexe Tsimembo Manambolomaty",
    nom == "Mandrozo" ~ "Zone Humide de Mandrozo",
    nom == "Paysage Harmonieux Protégés Bemanevika" ~ "Complexe des Zones Humides de Bemanevika",
    nom == "Nord Ifotaky" ~ "INord fotaky",
    TRUE ~ nom)) %>%
  arrange(nom_wdpa) %>%
  mutate(rownum = row_number())

# On sauvegarde le résultat
save(AP_Vahatra, file = "data/AP_Vahatra.rds")

# On ne garde que les aires de WDPA qui apparaissent dans Vahatra
WDPA_commun <- WDPA_Mada %>%
  filter(NAME %in% AP_Vahatra$nom_wdpa) %>%
  filter(!(NAME == "Analalava" & IUCN_CAT == "Not Reported")) %>%
  filter(!(NAME == "Site Bioculturel d'Antrema" & IUCN_CAT == "Not Reported")) %>%
  filter(DESIG != "UNESCO-MAB Biosphere Reserve") %>%
  arrange(NAME)  %>%
  mutate(rownum = row_number())
       
# Cette fonction calcule la part d'un polygone incluse dans un 
# autre polygone et retourne un ratio entre 0 et 1
ratio_inclus <- function(x, y) {
  inclus <- st_intersection(x, y)
  ratio <- st_area(inclus) / st_area(x)
  return(ratio)
}

# On calcule la part des polygones Vahatra incluse dans les polgones WDPA 
V_in_W <- map2_dbl(WDPA_commun$geometry, AP_Vahatra$geometry, ratio_inclus)
# Puis l'inverse
W_in_V <- map2_dbl(AP_Vahatra$geometry, WDPA_commun$geometry, ratio_inclus)
# On fait un facteur des deux
recoupement_mutuel <- V_in_W * W_in_V
# Qu'on ramène dans les jeux de données d'origine
WDPA_commun2 <- bind_cols(WDPA_commun, V_in_W = V_in_W, W_in_V = W_in_V,
                         recoupement_mutuel = recoupement_mutuel) %>%
  arrange(recoupement_mutuel, rownum)
AP_Vahatra2 <- bind_cols(AP_Vahatra, V_in_W = V_in_W, W_in_V = W_in_V,
                        recoupement_mutuel = recoupement_mutuel) %>%
  arrange(recoupement_mutuel, rownum)

# On prend maintenant les 5 les plus éloignés et on les visualise
min_recoup <- WDPA_commun2 %>%
  filter(row_number() <= 10) %>%
  select(nom_wdpa = NAME, rownum) %>%
  mutate(source = "WDPA") %>%
  bind_rows(select(filter(AP_Vahatra2, rownum %in% .$rownum), nom_wdpa, rownum)) %>%
  mutate(source = ifelse(is.na(source), "Vahatra", source))
tmap_mode("plot")
min_recoup %>%
  tm_shape() +
  tm_polygons() +
  tm_facets(by = c("nom_wdpa", "source")) +
  tm_layout(panel.label.size=3)

On peut également comparer ceux pour lesquels on a des différences de date ou de statut.

Code

# On garde seulement les métadonnées qu'on veut comparer
WDPA_a_comparer <- WDPA_commun %>% # On repart des AP communes
  st_drop_geometry() %>% # Plus besoin de spatial
  select(nom_wdpa = NAME, type_wdpa = INT_CRIT, cat_iucn_wdpa = IUCN_CAT,
         year_wdpa = STATUS_YR) # On ne garde que les colonnes à comparer

verif_meta_wdpa <-AP_Vahatra %>%
  st_drop_geometry() %>% # Pas besoin d'un jeu spatial
  select(nom:date_modification, nom_wdpa) %>% # colonnes à garder dans Vahatra
  # On renomme la catégorie IUCN de Vahatra et on code les NA comme dans WDPA
  mutate(cat_iucn = ifelse(is.na(cat__iucn), "Not Reported", cat__iucn)) %>%
  relocate(cat_iucn, .before = cat__iucn) %>% # Nouvelle colonne près de l'ancienne
  left_join(WDPA_a_comparer, by = "nom_wdpa") %>% # On rassemble Vahatra et WDPA
  select(-nom_wdpa, -cat__iucn) %>% # On enlève les colonnes inutiles
  # On compare les dates et statuts
  mutate(`Différence de date` = year(date_creation) != year_wdpa,
         `Différence de statut` = cat_iucn != cat_iucn_wdpa)

verif_meta_wdpa %>%
  summarise(`Nombre d'aires protégées comparées` = n(),
            `Différence de date` = sum(`Différence de date`),
            `Différence de statut` = sum(`Différence de statut`)) %>%
  gt() %>%
  tab_header(title = paste("Différences entre les données de WDPA et celles de",
                     "l'assciation Vahatra sur les aires protégées terrestres",
                     "à Madagascar"))

Différences entre les données de WDPA et celles de l'assciation Vahatra sur les aires protégées terrestres à Madagascar
Nombre d'aires protégées comparées	Différence de date	Différence de statut
98	68	40

Dans les cas qu’on peut comparer, les données de l’association Vahatra semblent plus fiables. On va donc privilégier l’utilisation de ces dernières.

On va également visualiser les aires de WDPA qui ne sont pas contenues dans Vahatra.

Code

WDPA_exclu <- WDPA_Mada %>%
  filter(!(NAME %in% AP_Vahatra$nom_wdpa))

tmap_mode("view")
WDPA_exclu %>%
  tm_shape() +
  tm_polygons(col = "IUCN_CAT")

Code

ratio_terrestre <- function(x) {
  inclus <- st_intersection(x, contour_mada$geometry)
  ratio <- st_area(inclus) / st_area(x)
  return(ratio)
}

# On crée un grand polygone avec toutes les AP dans Vahatra
AP_Vahatra_fusion <- st_union(AP_Vahatra)

# On calcule pour les aires protégées de WDPA qui ne sont pas dans Vahatra
# dans quelle mesure elles sont terrestres et pas superposées à d'autres AP
# déjà dans Vahatra
WDPA_exclu <- WDPA_exclu %>%
  filter(st_geometry_type(.) == "MULTIPOLYGON") %>%
  mutate(part_terrestre = map2_dbl(.$geometry, contour_mada$geometry, 
                                   ratio_inclus),
         part_deja_autre = map2_dbl(.$geometry, AP_Vahatra_fusion, 
                                    ratio_inclus))

# On garde celles qui sont au moins 25% terrestre et 75% pas superposées
WDPA_a_inclure <- WDPA_exclu %>%
  filter(part_terrestre >= 0.25 & part_deja_autre <= 0.25) %>%
  mutate(full_name = paste(INT_CRIT, NAME)) %>%
  select(nom = NAME, WDPAID, full_name, creation = STATUS_YR)

On a visiblement des aires protégées qu’il serait pertinent d’inclure et qui ne sont pas dans Vahatra.

1.4 Enjeux de fiabilité des données d’aires protégées

Important pour l’analyse : si périmètres pas juste => phénomènes de leakage, faux positifs ou faux négatifs.

Enjeu aussi des métadonnées : date ou type sont importants pour l’analyse et celle-ci perd en fiabilité si ces informations ne sont pas correctes.

Goodman, Steven M., Marie Jeanne Raherilalao, Sébastien Wohlhauser, Jean Clarck N. Rabenandrasana, Herivololona M. Rakotondratsimba, Fanja Andriamialisoa, and Malalarisoa Razafimpahanana. 2018. Les Aires Protégées Terrestres de Madagascar: Leur Histoire, Description Et Biote. Association Vahatra.